import numpy as np
import pandas as pd

def haversine(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    c = 2 * np.arcsin(np.sqrt(a))
    r = 6371
    return r * c


stations_data = pd.DataFrame({
    'Station_ID': [1, 2, 3, 4, 5, 6, 7, 8, 9],
    'Longitude': [110.125713, 110.08442, 110.029866, 109.962839, 109.956003, 109.920425, 109.839046, 109.823329, 109.767127],
    'Latitude': [32.815024, 32.771676, 32.748994, 32.743622, 32.812194, 32.856136, 32.860495, 32.847468, 32.807855]
})
 
demands_data = pd.DataFrame({
    'Demand_ID': range(1, 51),
    'Longitude': [
        110.1053385, 110.1147032, 110.0862574, 110.0435344, 110.0575508,
        110.0386243, 110.0115086, 110.0390602, 110.0246454, 110.0575847,
        109.9456331, 109.9612274, 109.94592, 109.9316682, 109.9245376,
        109.7087533, 109.7748005, 109.7475891, 109.7534532, 109.783015,
        109.7410728, 109.7554844, 109.7147417, 109.8807093, 109.8070677,
        109.9054481, 109.8954509, 109.8979229, 109.8942179, 109.8610985,
        109.8744682, 109.8338804, 109.870924, 109.8292467, 109.8711312,
        109.8813363, 109.978788, 109.8166563, 109.8151216, 109.885638,
        109.9890984, 109.9647812, 109.9303732, 109.9401099, 109.944496,
        109.979708, 109.976757, 109.94999, 109.973673, 109.967765
    ],
    'Latitude': [
        32.77881526, 32.75599834, 32.74905239, 32.74275416, 32.76712584,
        32.70855831, 32.72619993, 32.73965997, 32.72360718, 32.76553658,
        32.7526657, 32.72286471, 32.70899877, 32.73848444, 32.70740885,
        32.7815564, 32.80016336, 32.80903496, 32.85129032, 32.82296929,
        32.82914197, 32.80581363, 32.79995734, 32.89696579, 32.79622985,
        32.89437141, 32.86724756, 32.83444574, 32.83224374, 32.90687042,
        32.89939698, 32.85616627, 32.848223, 32.83825122, 32.88979101,
        32.8642824, 32.75943454, 32.8096699, 32.82822489, 32.84032485,
        32.80854774, 32.80993619, 32.78956582, 32.85264625, 32.802178,
        32.817449, 32.811064, 32.795207, 32.746858, 32.820998
    ],
    'Demand_kg': [3, 4, 2, 0, 8, 7, 4, 9, 10, 6, 7, 12, 3, 5, 6, 5, 3, 13, 12, 3,
                  14, 10, 4, 34, 6, 6, 3, 4, 20, 5, 6, 5, 3, 15, 2, 6, 3, 4, 3, 2,
                  6, 5, 9, 3, 3, 4, 6, 4, 4, 0]
})

stations_lons = stations_data['Longitude'].values
stations_lats = stations_data['Latitude'].values
demands_lons = demands_data['Longitude'].values
demands_lats = demands_data['Latitude'].values

stations_lons_matrix, demands_lons_matrix = np.meshgrid(stations_lons, demands_lons)
stations_lats_matrix, demands_lats_matrix = np.meshgrid(stations_lats, demands_lats)

distances_matrix = np.zeros((stations_lons_matrix.shape[0], demands_lons_matrix.shape[1]), dtype=float)

for i in range(stations_lons_matrix.shape[0]):
    for j in range(demands_lons_matrix.shape[1]):
        distance = haversine(
            stations_lons_matrix[i, j], stations_lats_matrix[i, j],
            demands_lons_matrix[i, j], demands_lats_matrix[i, j]
        )
        distances_matrix[i, j] = distance

print("距离矩阵 (所有点对之间的距离):")
print(distances_matrix)
